%
%		LAMA_Implementation_Doc
%
\documentclass[pdftex,a4paper,parskip,listof=totoc,bibliography=totoc,onehalfspacing,12pt]{scrreprt}
\usepackage[english]{babel} % If german is needed add ngerman in first bracket 
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{textcomp}
\usepackage{multirow}
%\usepackage{graphicx} is loaded with pdfpages
\usepackage{pict2e}
\usepackage{multicol}
\usepackage[automark]{scrpage2}
\usepackage{setspace}
\usepackage{color}
\usepackage{natbib}
\usepackage{booktabs}
\usepackage{amsmath}
\usepackage{mathtools}			%For use of \mathrlap and \mathllap
\usepackage{amsfonts}			%For use of \mathbb{R} (e.g. show real number letter)
\usepackage{float}
\usepackage{longtable}
\usepackage{pdfpages}
\usepackage{eso-pic}
\usepackage{dirtree}
\usepackage{multirow}
\usepackage{gauss}				% For use of \newcommand*\dashline... (vertical line in arrays)
\usepackage[vlined]{algorithm2e}
\usepackage{enumitem}
	\setlist[itemize]{itemsep=-12pt}
	\setlist[enumerate]{itemsep=-12pt}
\usepackage{subcaption}
\usepackage{lscape}		%For use of \begin{landscape} \end{landscape} Zur Darstellung einzelner Seiten im Querformat
\usepackage{adjustbox}
\usepackage{nicefrac}
\usepackage{physics}
\KOMAoptions{twoside}
\setcounter{MaxMatrixCols}{50}

\usepackage[left=25mm,right=25mm,bottom=35mm, footskip=20mm]{geometry}

\usepackage{siunitx}		% Show units
\sisetup{per-mode=fraction}	% Representation of fractions


\usepackage{listings}
\usepackage{scrhack}		% Comma-option \KOMAoption{listof}{leveldown} only effects lstlistoflistings 

\lstdefinestyle{MyMatlab}		% Definition of individual listing-styles
{						% Applications with "\lstinputlisting[style=MyMatlab]{../reflexion/Matlab_Skripte/resample.m}"
	language={Matlab},
	captionpos={t},
	inputencoding={latin1},
	showstringspaces={false},
	frame={tlRB},
	basicstyle=\scriptsize
}

\lstdefinestyle{MyC++}
{
	language={C++},
	tabsize={8},
	captionpos={t},
	frame={tlRB},
	basicstyle=\tiny
}

\lstdefinestyle{MyFortran}
{
	language={[77]Fortran},
	captionpos={t},
	frame={tlRB},
	basicstyle=\tiny
}

\lstdefinestyle{MySh}
{
	language={sh},
	captionpos={t},
	frame={tlRB},
	basicstyle=\tiny
}

\lstdefinestyle{MyGnuplot}
{
	language={Gnuplot},
	captionpos={t},
	frame={tlRB},
	basicstyle=\scriptsize
}

\lstset{literate=%
  {Ö}{{\"O}}1
  {Ä}{{\"A}}1
  {Ü}{{\"U}}1
  {ß}{{\ss}}2
  {ü}{{\"u}}1
  {ä}{{\"a}}1
  {ö}{{\"o}}1
}

\usepackage{rotating}
\usepackage{wrapfig}
\usepackage{tikz}
	\usetikzlibrary{shapes}
	\usetikzlibrary{decorations.pathmorphing}
	\usetikzlibrary{decorations.shapes}
	\usetikzlibrary{shapes.geometric}
	\pgfdeclarelayer{edgelayer}
	\pgfdeclarelayer{nodelayer}
	\pgfsetlayers{edgelayer,nodelayer,main}
	\tikzstyle{none}=[inner sep=0pt]
	\tikzset{decorate with/.style={decorate,decoration={shape backgrounds,shape=#1,shape size=2mm}}}

\definecolor{blue}{rgb}{0.01,0.01,0.95}	%HEX: 0202F2
\definecolor{green}{rgb}{0.01,0.65,0.01}	%HEX: 02A602
\definecolor{red}{rgb}{0.95,0.01,0.01}	%HEX: F20202
\definecolor{brown}{rgb}{0.55,0.27,0.07}	%HEX: 8B4513
\definecolor{violet}{rgb}{0.65,0.05,0.7}	%HEX: A60DB3

\setcounter{tocdepth}{1}	%Legt die Anzeigetiefe des Inhaltsverzeichnisses fest

\renewcommand{\theequation}{\arabic{section}.\arabic{equation}} %Die Nummerierung Formeln wird geändert


\renewcommand{\partheadmidvskip}{\enskip}
\renewcommand{\partformat}{\thepart\autodot}

\renewcommand{\arraystretch}{1.2} %Legt den Faktor fest, um den der Zeilenabstand innerhalb einer Tabelle oder eines Array gedeht wird
\newcommand*\dashline{\hspace{-0.3em}\rotatebox[origin=c]{90}{\scalebox{-0.3}[1]{$-~~-$}}\hspace{-0.3em}}
\newcommand*\longdashline{\hspace{-0.3em}\rotatebox[origin=c]{90}{\scalebox{-0.4}[1]{$-~~-~~-~~-$}}\hspace{-0.3em}}

%\setlength{\oddsidemargin}{20mm} 
%\setlength{\evensidemargin}{20mm} 

\renewcommand*{\chapterheadstartvskip}{\vspace*{-1\baselineskip}}	%legt Abstand vor chapter fest

\newcommand{\CC}[1][]{$\text{C\hspace{0ex}}^{_{_{_{++}}}}					%Write nicer ++ signs in C++
                      \ifthenelse{\equal{#1}{}}{}{\text{\hspace{-.625ex}#1}}$} 
\newcommand{\clang}[1][]{$\text{clang\hspace{0ex}}^{_{_{_{++}}}}		%Write nicer ++ signs in clang++
                      \ifthenelse{\equal{#1}{}}{}{\text{\hspace{-.25ex}(#1)}}$} 
\newcommand{\gCC}[1][]{$\text{g\hspace{0ex}}^{_{_{_{++}}}}				%Write nicer ++ signs in g++
                      \ifthenelse{\equal{#1}{}}{}{\text{\hspace{-.25ex}(#1)}}$} 
                      
\newcommand\coolover[2]{\mathrlap{\smash{%
\overbrace{\phantom{\begin{matrix} #2 %
\end{matrix}}}^{\mbox{$#1$}}}}#2}
                      
\newcommand{\shellcmd}[1]{\indent\indent\texttt{#1}}	%change font for shell entries  
\newcommand{\shellcmdline}[1]{\indent\indent\texttt{\quad#1}} 	%change font for command lines and indention


\makeatother
\usepackage{hyperref}
\definecolor{LinkColor}{rgb}{0,0,0.75}
\hypersetup{                    %Farben der Links im pdf werden festgelegt
colorlinks=true,
linkcolor=LinkColor,
citecolor=LinkColor,
filecolor=LinkColor,
menucolor=LinkColor,
pagecolor=LinkColor,
urlcolor=LinkColor}

\KOMAoption{listof}{leveldown} %Bewirkt, dass listoffigures, listoftables und lstlistoflistings nicht als chapter behandelt werden. Erschenungsbild im Inhaltsverzeichnis und im Dokument selbst (kein automatischer seitenumbruch am Ende der Verzeichnisse).

\pagestyle{scrheadings} %Legt die Art des Seitenformats für alle folgenden Seiten fest. scrheadings: Um bei zweiseitig formatierten Dokumenten (=Bücher, die unterschiedliche "`linke"' und "`rechte"' Seiten haben) den Seitenkopf und -fuß automatisch an die jeweilige Seite anzupassen, gibt es die folgenden Befehle: \ihead, \chead, \ohead, \ifoot, \cfoot, \ofoot für scrheadings wird "`\usepackage{scrpage2}"' benötigt
\clearscrheadfoot %Löscht alle Kopf- und Fußzeilen
\ohead{\headmark}
\automark[section]{chapter}
\ofoot[\pagemark]{\pagemark}		%[plain-Seiten]{normale Seiten}
\setheadsepline{0.4pt} %erzeugt eine Linie der Stärke 0.4pt zwischen Rumpf und Kopfzeile
\setlength{\headsep}{10mm} %Legt den Abstand zwischen der Kopfzeile und dem Rumpf der Seite fest

\hyphenation{CATIA} %definiert die Stellen zum Trennen von Wörter mit "`-"' Beispiel: \hyphenation{er-go-no-mic}

\title{LAMA_Implementation_Doc}
\author{}
\graphicspath{{@CMAKE_SOURCE_DIR@/../doc/guide/}}

% -----------------------------------------------------------------------------------------------
\begin{document}
% -----------------------------------------------------------------------------------------------

\selectlanguage{english}

\numberwithin{equation}{chapter} %legt die Formelnummerierung fest

\thispagestyle{empty} %Kopf und Fußzeile wird für diese Seite ausgeschaltet
\newgeometry{left=25mm,right=25mm,bottom=40mm,top=20mm}
\begin{figure}[h] % Gleiten ist durch Verwendung des H verhindert
\begin{flushright}
\includegraphics[scale=0.15]{./images/wave_logo.png}
\end{flushright}
\end{figure}

\begin{center}
\vspace{2cm}
\huge{WAVE-Inversion}\\
\vspace{0.5cm}
\large{HPC Full-Waveform Inversion Software}
\end{center}

\vfill
\begin{center}
{\Large{\url{wave-toolbox.org}}}


{\small Document created \today}
\end{center}

\newpage 
\thispagestyle{empty}
\begin{center}
{\large
Supported by the German Ministry of Education and Research (BMBF) through the project \textbf{WAVE}, grant 01IH15004A.
}
\end{center}

\cleardoublepage

\pagenumbering{Roman}
\setcounter{page}{1}
\restoregeometry


\newpage

\tableofcontents % Table of contents is inserted \tableofcontens (insert \addcontentsline, so the content is linked correctly)
\addcontentsline{toc}{chapter}{Contents} %"Contents" will be displayed in contents



\pagenumbering{arabic}
\setcounter{page}{1}

\cleardoublepage
%\markboth{}{License}
\chapter*{License}
\addcontentsline{toc}{chapter}{License}

WAVE-Inversion is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License (GNU AGPL v3) as published by the Free Software Foundation, version 3.0 of the License only.
 
WAVE-Inversion is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details. You should have received a copy of the GNU Affero General Public License along with FDSimulation. See file LICENSE and/or \url{https://www.gnu.org/licenses/agpl-3.0.de.html}.

The authors of WAVE-Inversion are listed in file \lstinline{AUTHORS}.

\chapter*{Acknowledgments}

We would like to thank the German Ministry of Education and Research (BMBF) who financially supported this full-waveform inversion toolbox through the project WAVE, grant 01IH15004A.

We would also like to acknowledge the developers of the FWI code IFOS2D (\url{http://www.gpi.kit.edu/Software-FWI.php}). IFOS2D served as a reference for our new implementation and we further used parts of the user manual of IFOS2D in this manual.

\chapter*{References}

A list of publications where WAVE-Inversion was used


% -----------------------------------------------------------------------------------------------
\cleardoublepage
\chapter*{Introduction}
\addcontentsline{toc}{chapter}{Introduction}
% -----------------------------------------------------------------------------------------------
 
The aim of full-waveform inversion (FWI) is to estimate the elastic material parameters in the underground. This can be achieved by minimizing the misfit energy between the modeled and field data using a gradient optimization
approach. Because FWI uses the full information content of each seismogram, structures below the seismic wave length can be resolved. This is a tremendous improvement in resolution compared to travel time tomography \citep{pratt:02}. The concept of FWI was originally developed by Albert Tarantola in the 1980s for the acoustic, isotropic elastic, and viscoelastic case \citep{tarantola:84a,tarantola:84b,tarantola:86,tarantola:88}. First numerical implementations were realized at the end of the 1980s \citep{gauthier:86,mora:87,pica:90}, but due to limited computational resources, the application was restricted to simple 2D synthetic test problems and small near offset datasets. At the begining of the 1990s the original time domain formulation was transfered to a robust frequency domain approach \citep{pratt:90a,pratt:90b}. With the increasing performance of supercomputers moderately sized problems could be inverted with frequency domain approaches. A spectacular result to prove the application of acoustic FWI on laboratory scale was presented by \cite{pratt:99} for ultrasonic tomography measurements on a simple block model. In a numerical blind test \cite{brenders:07} achieved a very good agreement between their inversion result and the unknown true P-wave velocity model. The parallelization and performance optimizations of the frequency domain approach \citep[see e.g. ][]{sourbier:09a,sourbier:09b} led to a wide range of acoustic FWI applications for problems on different scales, from the global scale, crustal scale over engineering and near surface scale, down to laboratory scale \citep{pratt:04}. Beside the application to geophysical problems, the acoustic FWI is also used to improve the resolution in medical cancer diagnostics \citep{pratt:07}. However, all these examples are restricted to the inversion of the acoustic material parameters: P-wave velocity, density and additionally the viscoacoustic damping $Q_p$ for the P-waves. Even today the independent 2D FWI of all three isotropic elastic material parameters is still a challenge. Most elastic approaches invert for P-wave velocity only and use empirical relationships to deduce the distribution of S-wave velocity and density \citep{shipp:02,sheen:06}. Recently some authors also investigated the independent multiparameter FWI in the frequency domain \citep{choi:08a,choi:08b,brossier:09}.
In order to extract information about the structure and composition of the crust from seismic observations, it is necessary to be able to predict how seismic wavefields are affected by complex structures. Since exact analytical
solutions to the wave equations do not exist for most subsurface configurations, the solutions can be obtained only by numerical methods. 

For iterative calculations of synthetic seismograms with limited computer resources fast and accurate modeling methods are needed. The FWI program WAVE-Inversion (Inversion of Full Observed Seismograms with C++) is using the finite difference simulation software WAVE-Simulation (insert link). It is also possible to use the WAVE-Inversion library with different simulation software.

% -----------------------------------------------------------------------------------------------
\cleardoublepage
\part{Theory of Full-Waveform Inversion}
\chapter{Derivation of the gradients of the objective function}
% -----------------------------------------------------------------------------------------------

% -----------------------------------------------------------------------------------------------
\section{Elastic}
\label{kap:el}
% -----------------------------------------------------------------------------------------------

3D elastic Wave Equation:
\begin{align}
\rho\pdv{v_x}{t}&=\pdv{\sigma_{xx}}{x}+\pdv{\sigma_{xy}}{y} + \pdv{\sigma_{xz}}{z}+ f_x\\
\rho\pdv{v_y}{t}&=\pdv{\sigma_{xy}}{x}+\pdv{\sigma_{yy}}{y} + \pdv{\sigma_{yz}}{z} + f_y\\
\rho\pdv{v_z}{t}&=\pdv{\sigma_{xz}}{x}+\pdv{\sigma_{yz}}{y} + \pdv{\sigma_{zz}}{z} + f_z\\
\pdv{\sigma_{xx}}{t} &= (\lambda+2\mu) \pdv{v_x}{x} + \lambda \pdv{v_y}{y} + \lambda \pdv{v_z}{z} + \pdv{\sigma_{xx0}}{t}\\
\pdv{\sigma_{yy}}{t} &= \lambda \pdv{v_x}{x} + (\lambda+2\mu) \pdv{v_y}{y} + \lambda \pdv{v_z}{z} + \pdv{\sigma_{yy0}}{t}\\
\pdv{\sigma_{zz}}{t} &= \lambda \pdv{v_x}{x} +  \lambda \pdv{v_y}{y} + (\lambda+2\mu) \pdv{v_z}{z} + \pdv{\sigma_{yy0}}{t}\\
\pdv{\sigma_{xy}}{t} &= \mu \left(\pdv{v_x}{y} +\pdv{v_y}{x}\right) + \pdv{\sigma_{xy0}}{t}\\
\pdv{\sigma_{yz}}{t} &= \mu \left(\pdv{v_y}{z} +\pdv{v_z}{y}\right) + \pdv{\sigma_{yz0}}{t}\\
\pdv{\sigma_{xz}}{t} &= \mu \left(\pdv{v_x}{z} +\pdv{v_z}{x}\right) + \pdv{\sigma_{xz0}}{t}
\end{align}

equation 4+5+6
\begin{equation}
 \pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}=\left(3\lambda+2\mu\right)\left(\pdv{v_x}{x}+\pdv{v_y}{y}+ \pdv{v_z}{z}\right)
\end{equation}


equation 4-6 rewritten:
\begin{align}
\pdv{v_x}{x}=\frac{1}{2\mu} \left(\pdv{\sigma_{xx}}{t} - \frac{\lambda}{3\lambda+2\mu} \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}\right)\right)\\
\pdv{v_y}{y}=\frac{1}{2\mu} \left(\pdv{\sigma_{yy}}{t} - \frac{\lambda}{3\lambda+2\mu} \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}\right)\right)\\
\pdv{v_z}{z}=\frac{1}{2\mu} \left(\pdv{\sigma_{zz}}{t} - \frac{\lambda}{3\lambda+2\mu} \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}\right)\right)
\end{align}

equation 4-6 rewritten:
\begin{align}
\pdv{v_x}{x}=A\pdv{\sigma_{xx}}{t} + B \pdv{\sigma_{yy}}{t}+B\pdv{\sigma_{zz}}{t}\\
\pdv{v_y}{y}=B\pdv{\sigma_{xx}}{t} + A \pdv{\sigma_{yy}}{t}+B\pdv{\sigma_{zz}}{t}\\
\pdv{v_z}{z}= B\pdv{\sigma_{xx}}{t} + B \pdv{\sigma_{yy}}{t}+A\pdv{\sigma_{zz}}{t}
\end{align}
With:
\begin{align}
A=\frac{1}{2\mu}-\frac{\lambda}{2\mu(N\lambda+2\mu)}\\
B=-\frac{\lambda}{2\mu(N\lambda+2\mu)}
\end{align}
With $N\in[2,3]$ (Number of dimensions)
Forward problem Matrix formulation;

\begin{equation}
 M^{-1}\left( \pdv{}{t}\mathbf{u}-b \right) =Q\mathbf{u}
\end{equation}


\begin{equation}
\mathbf{u}=\left(v_x,v_y,v_z,\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy},\sigma_{yz} ,\sigma_{xz} \right)^T
\end{equation}

\begin{equation}
\mathbf{b}=\left(f_x,f_y,f_z,\pdv{\sigma_{xx0}}{t},\pdv{\sigma_{yy0}}{t},\pdv{\sigma_{zz0}}{t},\pdv{\sigma_{xy0}}{t},\pdv{\sigma_{yz0}}{t},\pdv{\sigma_{xz0}}{t}\right)^T
\end{equation}

\begin{equation}
Q=
 \begin{pmatrix}
   0          & 0         & 0         & \pdv{}{x} & 0         & 0         & \pdv{}{y} & 0         & \pdv{}{z}\\
   0          & 0         & 0         & 0         & \pdv{}{y} & 0         & \pdv{}{x} & \pdv{}{z} & 0        \\
   0          & 0         & 0         & 0         & 0         & \pdv{}{z} & 0         & \pdv{}{y} & \pdv{}{x}\\
   \pdv{}{x}  & 0         & 0         & 0         & 0         & 0         & 0         & 0         & 0        \\
   0          & \pdv{}{y} & 0         & 0         & 0         & 0         & 0         & 0         & 0        \\
   0          & 0         & \pdv{}{z} & 0         & 0         & 0         & 0         & 0         & 0        \\
   \pdv{}{y}  & \pdv{}{x} &  0        & 0         & 0         & 0         & 0         & 0         & 0        \\
   0          & \pdv{}{z} & \pdv{}{y} & 0         & 0         & 0         & 0         & 0         & 0        \\
   \pdv{}{z}  & 0         & \pdv{}{x} & 0         & 0         & 0         & 0         & 0         & 0        \\
 \end{pmatrix}
\end{equation}

\begin{equation}
M=
 \begin{pmatrix}
   \rho       & 0         & 0         & 0         & 0         & 0         & 0         & 0         & 0       \\
   0          & \rho      & 0         & 0         & 0         & 0         & 0         & 0         & 0        \\
   0          & 0         & \rho      & 0         & 0         & 0         & 0         & 0         & 0        \\
   0          & 0         & 0         & A         & B         & B         & 0         & 0         & 0        \\
   0          & 0         & 0         & B         & A         & B         & 0         & 0         & 0        \\
   0          & 0         & 0         & B         & B         & A         & 0         & 0         & 0        \\
   0          & 0         & 0         & 0         & 0         & 0         & 1/\mu      & 0         & 0        \\
   0          & 0         & 0         & 0         & 0         & 0         & 0         & 1/\mu      & 0        \\
   0          & 0         & 0         & 0         & 0         & 0         & 0         & 0         & 1/\mu     \\
 \end{pmatrix}
\end{equation}

\vspace{0.3cm}

\begin{center} \Large{Gradient:}\end{center}

\begin{equation}
 \pdv{E}{m}=\int_T dt\mathbf{\Psi}^T \pdv{M^{-1}}{m}\left( \pdv{}{t}\mathbf{u}-b \right)
\end{equation}

\begin{equation}
\mathbf{ \Psi}=\left(\psi_x,\psi_y,\psi_z,\phi_{xx},\phi_{yy},\phi_{zz},\phi_{xy},\phi_{yz} ,\phi_{xz} \right)^T
\end{equation}

\begin{align}
 \pdv{E}{\rho}&=\int_T dt \left[ \left(\pdv{v_x}{t}-f_x\right)\psi_x +\left(\pdv{v_y}{t}-f_y\right)\psi_y +\left(\pdv{v_z}{t}-f_z\right)\psi_z \right]\\
 \pdv{E}{\lambda}&=\pdv{A}{\lambda}\int_T dt  \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}\right) \left(\phi_{xx}+\phi_{yy}+\phi_{zz}\right) \\
 \pdv{E}{\mu}&=\int_T dt \left[\pdv{A}{\mu}\left(\pdv{\sigma_{xx}}{t}\phi_{xx}+\pdv{\sigma_{yy}}{t}\phi_{yy}+\pdv{\sigma_{zz}}{t}\phi_{zz}\right)\right. \\
 &+ \pdv{B}{\mu} \left(\left(\pdv{\sigma_{yy}}{t}+\pdv{\sigma_{zz}}{t}\right)\phi_{xx}
 + \pdv{B}{\mu} \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{zz}}{t}\right)\phi_{yy}
  + \pdv{B}{\mu} \left(\pdv{\sigma_{xx}}{t}+\pdv{\sigma_{yy}}{t}\right)\phi_{zz}\right)\\
  &-\left.\frac{1}{\mu^2}\left(\pdv{\sigma_{xy}}{t}\phi_{xy}+\pdv{\sigma_{yz}}{t}\phi_{yz}+\pdv{\sigma_{xz}}{t}\phi_{xz}\right)\right]
\end{align}

With:
\begin{align}
 \pdv{A}{\lambda}&=\pdv{B}{\lambda}=-\frac{1}{(N\lambda+2\mu)^2} \\
 \pdv{B}{\mu}&=\frac{N\lambda^2+4\mu\lambda}{2\mu^2(N\lambda+2\mu)^2} \\
 \pdv{A}{\mu}&=-\frac{1}{2\mu^2}+\pdv{B}{\mu}
\end{align}


% -----------------------------------------------------------------------------------------------
\section{Acoustic}
\label{kap:ac}
% -----------------------------------------------------------------------------------------------

3D acoustic Wave Equation:

\begin{align}
\rho\pdv{v_x}{t}&=\pdv{p}{x}+ f_x\\
\rho\pdv{v_y}{t}&=\pdv{p}{y} + f_y\\
\rho\pdv{v_z}{t}&=\pdv{p}{z} + f_z\\
\frac{1}{\lambda}\pdv{p}{t} &= \pdv{v_x}{x} +  \pdv{v_y}{y} + \pdv{v_z}{z} + \pdv{p}{t}\\
\end{align}

\begin{equation}
 M^{-1}\left( \pdv{}{t}\mathbf{u}-b \right) =Q\mathbf{u}
\end{equation}

\begin{equation}
\mathbf{u}=\left(v_x,v_y,v_z,p \right)^T
\end{equation}

\begin{equation}
\mathbf{b}=\left(f_x,f_y,f_z,\pdv{p}{t})^T\right)
\end{equation}

\begin{equation}
Q=
 \begin{pmatrix}
   0         & 0         & 0         & \pdv{}{x}\\
   0         & 0         & 0         & \pdv{}{y}\\
   0         & 0         & 0         & \pdv{}{z}\\
   \pdv{}{x} & \pdv{}{y} & \pdv{}{z} & 0        \\

 \end{pmatrix}
\end{equation}

\begin{equation}
M=
 \begin{pmatrix}
   \rho       & 0         & 0         & 0                \\
   0          & \rho      & 0         & 0                \\
   0          & 0         & \rho      & 0                \\
   0          & 0         & 0         & \frac{1}{\lambda}\\
 \end{pmatrix}
\end{equation}

\vspace{0.3cm}

\begin{center} \Large{Gradient:}\end{center}

\begin{equation}
 \pdv{E}{m}=\int_T dt\mathbf{\Psi}^T \pdv{M^{-1}}{m}\left( \pdv{}{t}\mathbf{u}-b \right)
\end{equation}

\begin{equation}
\mathbf{ \Psi}=\left(\psi_x,\psi_y,\psi_z,\phi\right)^T
\end{equation}

\begin{align}
 \pdv{E}{\rho}&= (-)\int_T dt \left[ \left(\pdv{v_x}{t}-f_x\right)\psi_x +\left(\pdv{v_y}{t}-f_y\right)\psi_y +\left(\pdv{v_z}{t}-f_z\right)\psi_z \right]\\
 \pdv{E}{\lambda}&=- \frac{1}{\lambda^2}\int_T dt  \pdv{p}{t}\phi \\
\end{align}



\part{Software Guide to WAVE-Inversion}
\chapter{Software Requirements and Quick Guide}

WAVE-Inversion is developed as a part of the WAVE-project (\url{http://wave-toolbox.org}) and is a seismic full-waveform inversion framework. 

\section{Software Requirements}

To begin with, a \CC-compiler has to be available. A compiler with  \CC[11]-standard support is needed. To compile WAVE-Inversion, make sure a couple of libraries are installed and running properly. All necessary libraries for proper usage and a short description can be seen in table \ref{tab:pack}. A \CC-compiler and LAMA are mandatory whereas the other libraries and software is optional to execute WAVE-Inversion.
As the software is meant for high performance computing (HPC) on cluster networks, it is primarily designed to run on Linux-platforms. The code has been tested using \gCC[4.9.3], \clang[3.9] as a  \CC-compiler and using OpenMPI (1.7.2, 1.10.0) as MPI implementation.

\begin{table}[h!]
\caption{List of necessary libraries, packages and compiler.}\label{tab:pack}
\centering
\begin{adjustbox}{width=1\textwidth,center=\textwidth}
	\begin{tabular}{lll}
	\toprule
	Program & Brief description & Weblink \\
	\midrule
	\CC-compiler	& The code is written in \CC[], the compiler	has 		& \\
										& to be compatible to compile \CC[11]-standard					& \\
LAMA 				& The LAMA-framework for high performance			& \url{https://www.libama.org} \\
	& computing 	(HPC)							    & \\
	\midrule
	OpenMPI			& Implementation of MPI for parallelised  			& \url{https://www.open-mpi.org} \\
										& calculations						& \\
		boost					& Implementation of the boost-library	for  				& \url{https://www.boost.org}		\\
										& debugging																& \\										
	Python				& Models can be generated using Python 		& \url{https://www.python.org}\\
	Matlab				& Models can be generated using Matlab		& \url{http://www.mathworks.com}\\	
	
	
	Seismic Un*x	& 	Seismic processing package, to generate a  	& \url{http://www.cwp.mines.edu/cwpcodes}\\ 
											&  Seismogram from SU-data  															& \\
	
	
	Sphinx				& 	A Python documentation generator 							& \url{http://www.sphinx-doc.org}\\
	doxygen			& 	Generation of documentation from annotated				& \url{http://www.doxygen.org}\\
										& \CC-code 															& \\
	googletest & Implementation of gtest for unit-tests & \url{http://code.google.com/p/googletest}\\
	\bottomrule
	\end{tabular}
	\end{adjustbox}
\end{table}

\section{Quick Guide}

\subsection{Installation}
To get started Download, WAVE-Inversion. 

You can either download the software directly from \url{https://git.scc.kit.edu/WAVE/WAVE-Inversion/repository/archive.zip?ref=master} and unzip the folder (e.g. by entering \shellcmd{tar -zxvf WAVE-Inversion.tgz}).

Alternatively you can also get WAVE-Inversion by entering the command: \\\shellcmdline{git clone https://git.scc.kit.edu/WAVE/WAVE-Inversion.git}

As preparatory steps environment variables have to be linked to LAMA. Add the following lines to the \shellcmd{{\char`\~}/.bashrc}:
\\\shellcmdline{ `export SCAI\_ROOT=[PATH\_TO\_LAMA\_BUILD]`}
\\\shellcmdline{ `export DYLD\_LIBRARY\_PATH=\$\{SCAI\_ROOT\}/lib:\$\{DYLD\_LIBRARY\_PATH\}`}
\\\shellcmdline{ `export LD\_LIBRARY\_PATH=\$\{SCAI\_ROOT\}/lib:\$\{LD\_LIBRARY\_PATH\}`}
\\\shellcmdline{ `export OMP\_NUM\_THREADS=1`}\\
You can also find the instructions in \shellcmd{WAVE-Inversion/par/README.md}. 

Change to the directory: \\\shellcmdline{build/}

You can compile the source code by simply entering: \\\shellcmdline{make prog}\\ in your command line. \gCC is set to be the default \CC-compiler. To use an alternative compiler, the compiler setting can be change in the \shellcmd{Makefile}, or otherwise by executing:\\\shellcmdline{make CXX=<yourCompiler> prog}

The runtime configuration (such as the FD grid) can be changed in the \shellcmd{configuration.txt} file, located in the directory \shellcmd{par/configuration/}. A list of all input variables and their dependencies can be found in section \ref{sec:config}. 

To create a html and LaTeX or Pdf documentation of the software code itself, you can generate a doxygen documentation. Therefore change the directory to \shellcmd{build/} and run the command \shellcmd{make doc}. After running doxygen, the html, LaTeX and Pdf version of the reference manual can be found in \shellcmd{build/doc/doxygen/} and \shellcmd{build/doc/guide/}.

\subsection{Example}

To run the simulation change the directory to: \\\shellcmdline{par/}\\
and start the program by running the shell script \shellcmd{start\_Inversion.sh} with  \\\shellcmdline{source start\_Inversion.sh}\\

The synthetic Seismogram data is saved in \shellcmd{par/seismograms/} and can be plotted by using either the \shellcmd{plotSeismogram.m} Matlab-script or by executing the Python-script \shellcmd{plotSeismogram.py}.
As a first test, run the default example by following the commands to obtain the seismogram. You can compare the seismogram to the seismogram shown in figure (add figure).


\section{Configuration}\label{sec:config}

\begin{table}[h!]
\caption[List of all configuration parameters.]{List of all configuration parameters, that can be added and changed in the config-file.}\label{tab:config}
\centering
\begin{adjustbox}{width=1\textwidth,center=\textwidth}
	\begin{tabular}{llll}
	\toprule
         Variable                 & Short description                                                   & Type   & Example value \\
	\midrule
         fieldSeisName            & Name of field data                                                  & string & ci/rectangle.true  \\    
         writeGradient            & Write gradient to file                                              &  int   & 1 (=yes) \\
         writeGradientPerShot     & Write gradient of each shot separately to file                      &  int   & 0 (=no) \\
         gradientFilename         & Name of gradients                                                   & string & gradients/grad \\     
         logFilename              & Name of log file                                                    & string & logs/steplengthSearch.log  \\
         verbose                  & only in lbfgs branch!, use detailed output                          &  int   & 0 (=no) \\
        \midrule
         maxIterations            & Maximum number of inversion iterations (per workflow stage if used) &  int   & 20  \\                 
         misfitType               & Type of misfit                                                      & string & L2  \\
         optimizationType         & Type of optimization                                                & string & conjugateGradient \\
         workflowFilename         & Name of workflow file                                               & string & workflow/workflow.txt \\
        \midrule
         steplengthInit           & Initial step length used in the first iteration of each stage       & double & 0.03  \\
         steplengthMin            & Minimum step length                                                 & double & 0.001 \\
         steplengthMax            & Maximum step length                                                 & double & 0.1 \\                     
         maxStepCalc              & Maximum number of additional step length calculations               &  int   & 4 \\                         
         scalingFactor            & Factor for multiplication or division of test step length           & double & 2.0 \\                       
         testShotStart            & Shot number of first test shot                                      &  int   & 0 \\                        
         testShotEnd              & Shot number of last test shot                                       &  int   & 17 \\                        
         testShotIncr             & Increment of test shots                                             &  int   & 1 \\                       
        \midrule
         useSourceSignalInversion & Use source time inversion                                           &  int   & 0 (=no) \\
         waterLevel               & Water level of source time inversion                                & double & 0.01 \\                        
        \midrule
         sourceTaperRadius        & Circular source taper: Radius in grid points                        &  int   & 20 \\ 
         receiverTaperRadius      & Circular receiver taper: Radius in grid points                      &  int   & 20 \\
         useEnergyPreconditioning & Approximated diagonal of Hessian is applied to gradient per shot    &  int   & 1 (=yes) \\ 
         epsilonHessian           & Water level to stabilize matrix inversion (recommended: 0.005)      & double & 0.005 \\      
         saveApproxHessian        & Save approximated diagonal of Hessian                               &  int   & 0 (=no) \\
         approxHessianName        & Name of file for approximated diagonal of Hessian                   & string & gradients/Hessian \\
	\bottomrule
	\end{tabular}
	\end{adjustbox}
\end{table}	


As mentioned in the quick guide, the configuration can be found in \shellcmd{par/configuration/}. In the \shellcmd{configuration.txt} file Input variables can be added or changed. A list of all variables and config parameters that can be added and changed can be found in table \ref{tab:config}. In the following subsection you can find a detailed description of all parameters.

\subsection{General in- and output}

Full-waveform inversion is an algorithm which needs and produces a lot of data. In this part, the general in- and output is described. Some in-/output options are very specific and therefore described in the corresponding section.

The core of every FWI is the field data that should be inverted. The folder and basename of the field data has to be defined with the parameter \verb+fieldSeisName+. The program is automatically appending the shot number and the seismogram type (p,vx,vy,vz). For example, if the parameter \verb+fieldSeisName+ is set to \verb+seismograms/fieldData+ and pressure seismograms are inverted, the complete file name that is tried to be read is \verb+seismograms/fieldData.shot_0.p.mtx+. Shots are counted from 0.
Synthetic seismograms are saved in the same way with the filename \verb+SeismogramFilename+. Note that seismograms with \verb+.It_0.+ correspond to the seismograms of the initial model of each workflow stage.

With the parameters \verb+writeGradientPerShot+ and \verb+writeGradient+, it can be decided whether the gradients per shot and the summed gradients over all shots should be saved to disc. If yes, the gradients are stored with the name \verb+gradientFilename+. The workflow stage and iteration number and the parameter class are automatically appended, e.g. with \verb+gradients/grad+ as \verb+gradientFilename+ the first gradient with respect to the P-wave velocity will be \verb+gradients/grad.stage_1.It_1.vp.mtx+.

The program will create a log file with the name specified by \verb+logFilename+. The file contains the misfit evolution (last column) for each stage and iteration (first two columns). The third column (optimum step length) is the step length which was used for the model update. The next three columns correspond to the three step length of the inexact line search and the subsequent three columns are the corresponding misfits. Note that these three misfits could be calculated only for a subset of the shots depending on the setting. Rows with iteration 0 denote the initial misfit of the corresponding workflow stage.

Additional to the configuration file a workflow file has to be used because some parameters are only defined there and are not available in the configuration file! The name is set by the parameter \verb+workflowFilename+. At the moment, the workflow file contains the following parameters: \verb+invertForVp+, \verb+invertForVs+, \verb+invertForDensity+, \verb+relativeMisfitChange+, \verb+filterOrder+, \verb+lowerCornerFreq+, \verb+upperCornerFreq+. With the first three parameters it can be decided which model parameter class should be updated. The parameter \verb+relativeMisfitChange+ will be explained in subsection \ref{config:abort} and the next three parameters are related to the frequency filtering which is documented in the manual of WAVE-Simulation (insert link). 

The \verb+verbose+ option for detailed output is not available yet, but will be implemented in the future.


\subsection{Misfit definition}

The misfit definition has to be specified with the parameter \verb+misfitType+. Currently, only L2 misfit is available and can be chosen by setting \verb+misfitType+ to \verb+L2+. Note that seismograms can be normalized for the calculation of the misfit and the adjoint sources. This option is recommended for field data.

\subsection{Step length estimation}

For gradient-based, local optimization algorithms usually a step length estimation is necessary. In the case the steepest descent or conjugate gradient method is used (with or without preconditioning), WAVE-Inversion performs an inexact line search using a parabolic fit of the misfit function. For the parabolic fit three step length/misfit pairs are used. To save computation time the first pair is given by the step length zero and the misfit of the current model. The second pair is given by the parameter \verb+steplengthInit+ and the misfit calculated from the updated model with \verb+steplengthInit+. For the third pair two cases are distinguished: 1) if the misfit was decreased, the step length is multiplied with the parameter \verb+scalingFactor+ and a third step length is searched which gives a larger misfit than the second step length. If the misfit is still decreasing, the step length will again be multiplied with \verb+scalingFactor+ until it increases or the parameter \verb+maxStepCalc+ is reached. 2) If the misfit was increased, the step length is divided by the parameter \verb+scalingFactor+ and a third step length is searched which gives a smaller misfit than the first step length. If the misfit is too large, the step length will again be divided by \verb+scalingFactor+ until it sufficiently decreases or the parameter \verb+maxStepCalc+ is reached.
Note that the parameter \verb+steplengthInit+ is decreased by 2\% after each iteration within each workflow stage. The parameter is reset at the beginning of each stage to the value specified in the configuration file.
After calculating three step length/misfit pairs, four cases are distinguished to choose a proper step length:

\begin{enumerate}
 \item $E_2$ < $E_1 \land E_3 > E_2$: apply parabolic fit to find step length or use \verb+steplengthMax+ 
 \item $E_2$ < $E_1 \land E_3 < E_2$: use third step length or \verb+steplengthMax+
 \item $E_2$ > $E_1 \land E_3 < E_1$: use third step length 
 \item $E_2$ > $E_1 \land E_3 > E_1$: use \verb+steplengthMin+
\end{enumerate}

Here, $E$ stands for misfit. Note that the parameter \verb+steplengthMax+ is used if the third step length (cases 1 and 2) exceeds the specified value.
To save computation time, the step length estimation can be performed with a subset of the sources. The sources that are used can be specified with the parameters \verb+testShotStart+, \verb+testShotIncr+ and \verb+testShotEnd+. The first shot that is used is \verb+testShotStart+. It is incremented by \verb+testShotIncr+ until it reaches \verb+testShotEnd+.

Currently it is only possible to use one common step length for all model parameter classes (P-wave velocity, S-wave velocity, etc.). In the case of steepest descent or conjugate gradient method (with or without preconditioning) the model is updated in the following way
\begin{equation}
 \vec{m}_{k+1} = \vec{m}_k - \alpha \cdot \frac{\mathrm{max}(\vec{m}_k)}{\mathrm{max}(\nabla_{\vec{m}} E_k)} \nabla_{\vec{m}} E_k   
\end{equation}
with $k$ as iteration number and $\alpha$ as step length. Thus, a step length of 0.01 results in a maximum model update of 1\% of the maximum of the corresponding model parameter.

\subsection{Gradient preconditioning}
\label{config:precond}

Usually it is helpful to manipulate the gradient, e.g. to mitigate artifcats close to the source and receiver positions. There are currently three possibilities for gradient preconditioning: source and receiver tapers and energy preconditioning.

Circular tapers around the source and receiver positions can be used by setting \verb+sourceTaperRadius+ > 0 and \verb+receiverTaperRadius+ > 0, respectively. The number corresponds to the radius of the taper in grid points. The taper is zero in the middle and increases logarithmically to one at the outside of the circle. 

Energy preconditioning can be used by setting \verb+useEnergyPreconditioning+ to 1. It is based on migration weight $K^{(1)}$ from \cite{plessix:04} (see also \cite{shin:01}) which is the inverse of an approximation of the diagonal of the Hessian. This type of preconditioning has also the effect of mitigating source/receiver artifacts. Moreover, it increases model updates in less illuminated areas of the model and can improve the convergence behaviour of the inversion. It is calculated by squaring and adding all available velocity wavefields summed over all time steps. The taper is calculated and applied for each shot separately. To stabilize the approximated Hessian a water level has to be set with the parameter \verb+epsilonHessian+. From experience, we recommend a value of 0.005.
The taper can be saved to disc with the file name \verb+approxHessianName+ by setting \verb+saveApproxHessian+ to 1.


\subsection{Optimization}

Currently, there are two different optimization methods available, the steepest descent and the conjugate gradient method. Both can be used with gradient preconditioning (see subsection \ref{config:precond}). The method has to be chosen with the parameter \verb+optimizationType+. Possible values are \verb+steepestDescent+ and \verb+conjugateGradient+. The characters are internally transformed to lowercase letters, so \verb+STEEPESTDESCENT+ would also be valid.

The conjugate gradient direction at iteration $k$, $c_k$, is calculated in the following way:
\begin{equation*}
 c_k =  \nabla_{\vec{m}} E_k + \beta_k \cdot c_{k-1}
\end{equation*}

The weight, $\beta_k$, is calculated after Polak and Ribi\'{e}re:
\begin{equation*}
 \beta_k^{\mathrm{PR}} = \frac{ \nabla_{\vec{m}}^T E_k (\nabla_{\vec{m}} E_k - \nabla_{\vec{m}} E_{k-1}) }{\nabla_{\vec{m}}^T E_{k-1} \nabla_{\vec{m}} E_{k-1}}
\end{equation*}

Note that a steepest descent update is performed in the first iteration of every workflow stage.

\subsection{Source time function inversion}

\subsection{Abort criterion}
\label{config:abort}

Currently, there are two abort criteria implemented. The first criterion determines the maximum iteration number per workflow stage which can be set with the parameter \verb+maxIterations+. The second criterion measures the relative misfit change, $\Delta E$, between the current iteration and the second last iteration:

\begin{equation*}
 \Delta E = \frac{|E_k-E_{k-2}|}{E_{k-2}}
\end{equation*}

If this value is smaller than the parameter \verb+relativeMisfitChange+ the workflow stage is changed or - if it is the last stage - the inversion stops.


\clearpage
\section{Pre- and Post-Processing}\label{sec:process}

As already mentioned in section \ref{sec:config}, data can be read from a single file as well as from a partitioned file-block. The advantage of a distributed file system ist, that CPUs/GPUs can read data simultaneously and therefore with substantial timesaving. Installing LAMA with all examples (\shellcmd{-DBUILD\_EXAMPLES=ON}), LAMA includes examples to partition and repartition files or file-blocks in \shellcmd{.mtx}-format (\shellcmd{vectorRepartition.exe}). 

To partition a single vector-file change to the directory where the files are located and run the following comand:  \\
\shellcmdline{\$SCAI\_ROOT/lama/examples/io/vectorRepartition.exe filename.mtx 1 }\\
\shellcmdline{filename\_\%r.mtx <NProcessors>}\\
The first part marks the LAMA-example to repartition a file, \shellcmd{filename.mtx} is a single file (\shellcmd{1}) you want to partition and \shellcmd{filename\_\%r.mtx} is the name of the output file-block. The number of files in a file-block has to be the same as the number of CPUs/GPUs (\shellcmd{<NProcessors>}) reading the file-block later on. 

To merge a fileblock \shellcmd{filename\_\%r.mtx} to a single vector-file, the command is very similar:\\
\shellcmdline{\$SCAI\_ROOT/lama/examples/io/vectorRepartition.exe filename\_\%r.mtx}\\
\shellcmdline{<NProcessors> filename.mtx 1}\\
In this case the name of the file-block and its size is written fist and then the name of the single output-file.

\section{Unit Tests}
To guarantee a functioning framework, WAVE-Inversion includes unit-tests using the googletest framework. To run all unit-tests, you first have to install googletest (e.g. from \url{http://code.google.com/p/googletest}). 

As preparatory steps environment variables have to be linked to googletest. Add the following lines to the \shellcmd{{\char`\~}/.bashrc}:
\\\shellcmdline{ `export GTEST\_DIR=[PATH\_TO\_GTEST\_DIRECTORY]`}
\\\shellcmdline{ `export LD\_LIBRARY\_PATH=\$\{GTEST\_DIR\}:\$\{LD\_LIBRARY\_PATH\}`}
\\\shellcmdline{ `export DYLD\_LIBRARY\_PATH=\$\{GTEST\_DIR\}:\$\{DYLD\_LIBRARY\_PATH\}}

To compile all tests, change the directory: \\\shellcmdline{/src/} \\ and compile the tests by simply entering:\\\shellcmdline{make utest} and \shellcmdline{make itest}

The executable file is located in the directory:
\\\shellcmdline{/build/bin}\\
To execute the unit-tests change to the \shellcmd{/par} folder and run the \shellcmd{Test\_unit} file by entering:
\\\shellcmdline{./../bin/Tests/Test\_unit}\\ 
For faultless function of WAVE-Inversion all tests should execute without failure.

\section{Benchmark}

\cleardoublepage
\addtocontents{toc}{\protect\setcounter{tocdepth}{0}}
\listoffigures
\addtocontents{toc}{\protect\setcounter{tocdepth}{1}}
\addcontentsline{toc}{chapter}{List of Figures}
\addtocontents{toc}{\protect\setcounter{tocdepth}{0}}
\listoftables
\addtocontents{toc}{\protect\setcounter{tocdepth}{1}}
\addcontentsline{toc}{chapter}{List of Tables}
\cleardoublepage
\bibliography{@CMAKE_SOURCE_DIR@/../doc/guide/WAVE_guide}
\bibliographystyle{apalike}
\end{document}
